home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DMTDS.FOR < prev    next >
Text File  |  1985-11-29  |  5KB  |  157 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DMTDS
  5. C
  6. C        PURPOSE
  7. C           MULTIPLY A GENERAL MATRIX A ON THE LEFT OR RIGHT BY
  8. C           INVERSE(T),INVERSE(TRANSPOSE(T)) OR INVERSE(TRANSPOSE(T*T))
  9. C           THE TRIANGULAR MATRIX T IS STORED COLUMNWISE IN COMPRESSED
  10. C           FORM, I.E. UPPER TRIANGULAR PART ONLY.
  11. C
  12. C        USAGE
  13. C           CALL DMTDS(A,M,N,T,IOP,IER)
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           A     - GIVEN GENERAL MATRIX WITH  M ROWS AND N COLUMNS.
  17. C                   A MUST BE OF DOUBLE PRECISION
  18. C           M     - NUMBER OF ROWS OF MATRIX A
  19. C           N     - NUMBER OF COLUMNS OF MATRIX A
  20. C           T     - GIVEN TRIANGULAR MATRIX STORED COLUMNWISE UPPER
  21. C                   TRIANGULAR PART ONLY. ITS NUMBER OF ROWS AND
  22. C                   COLUMNS K IS IMPLIED BY COMPATIBILITY.
  23. C                   K = M IF IOP IS POSITIVE,
  24. C                   K = N IF IOP IS NEGATIVE.
  25. C                   T OCCUPIES K*(K+1)/2 STORAGE POSITIONS.
  26. C                   T MUST BE OF DOUBLE PRECISION
  27. C           IOP   - INPUT VARIABLE FOR SELECTION OF OPERATION
  28. C                   IOP = 1 - A IS REPLACED BY INVERSE(T)*A
  29. C                   IOP =-1 - A IS REPLACED BY A*INVERSE(T)
  30. C                   IOP = 2 - A IS REPLACED BY INVERSE(TRANSPOSE(T))*A
  31. C                   IOP =-2 - A IS REPLACED BY A*INVERSE(TRANSPOSE(T))
  32. C                   IOP = 3 - A IS REPLACED BY INVERSE(TRANSPOSE(T)*T)*A
  33. C                   IOP =-3 - A IS REPLACED BY A*INVERSE(TRANSPOSE(T)*T)
  34. C           IER   - RESULTING ERROR PARAMETER
  35. C                   IER =-1 MEANS M AND N ARE NOT BOTH POSITIVE
  36. C                                 AND/OR IOP IS ILLEGAL
  37. C                   IER = 0 MEANS OPERATION WAS SUCCESSFUL
  38. C                   IER = 1 MEANS TRIANGULAR MATRIX T IS SINGULAR
  39. C
  40. C        REMARKS
  41. C           SUBROUTINE DMTDS MAY BE USED TO CALCULATE THE SOLUTION OF
  42. C           A SYSTEM OF EQUATIONS WITH SYMMETRIC POSITIVE DEFINITE
  43. C           COEFFICIENT MATRIX. THE FIRST STEP TOWARDS THE SOLUTION
  44. C           IS TRIANGULAR FACTORIZATION BY MEANS OF DMFSD, THE SECOND
  45. C           STEP IS APPLICATION OF DMTDS.
  46. C           SUBROUTINES DMFSD AND DMTDS MAY BE USED IN ORDER TO
  47. C           CACULATE THE PRODUCT TRANSPOSE(A)*INVERSE(B)*A WITH GIVEN
  48. C           SYMMETRIC POSITIVE DEFINITE B AND GIVEN A IN THREE STEPS
  49. C           1) TRIANGULAR FACTORIZATION OF B (B=TRANSPOSE(T)*T)
  50. C           2) MULTIPLICATION OF A ON THE LEFT BY INVERSE(TRANSPOSE(T))
  51. C              A IS REPLACED BY C=INVERSE(TRANSPOSE(T))*A
  52. C           3) CALCULATION OF THE RESULT FORMING TRANSPOSE(C)*C
  53. C
  54. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  55. C           NONE
  56. C
  57. C        METHOD
  58. C           CALCULATION OF X = INVERSE(T)*A IS DONE USING BACKWARD
  59. C           SUBSTITUTION TO OBTAIN X FROM T*X = A.
  60. C           CALCULATION OF Y = INVERSE(TRANSPOSE(T))*A IS DONE USING
  61. C           FORWARD SUBSTITUTION TO OBTAIN Y FROM TRANSPOSE(T)*Y = A.
  62. C           CALCULATION OF Z = INVERSE(TRANSPOSE(T)*T)*A IS DONE
  63. C           SOLVING FIRST TRANSPOSE(T)*Y = A AND THEN T*Z = Y, IE.
  64. C           USING THE ABOVE TWO STEPS IN REVERSE ORDER
  65. C
  66. C     ..................................................................
  67. C
  68.       SUBROUTINE DMTDS(A,M,N,T,IOP,IER)
  69. C
  70. C
  71.       DIMENSION A(1),T(1)
  72.       DOUBLE PRECISION DSUM,A,T
  73. C
  74. C        TEST OF DIMENSION
  75.       IF(M)2,2,1
  76.     1 IF(N)2,2,4
  77. C
  78. C        ERROR RETURN IN CASE OF ILLEGAL DIMENSIONS
  79.     2 IER=-1
  80.       RETURN
  81. C
  82. C        ERROR RETURN IN CASE OF SINGULAR MATRIX T
  83.     3 IER=1
  84.       RETURN
  85. C
  86. C        INITIALIZE DIVISION PROCESS
  87.     4 MN=M*N
  88.       MM=M*(M+1)/2
  89.       MM1=M-1
  90.       IER=0
  91.       ICS=M
  92.       IRS=1
  93.       IMEND=M
  94. C
  95. C        TEST SPECIFIED OPERATION
  96.       IF(IOP)5,2,6
  97.     5 MM=N*(N+1)/2
  98.       MM1=N-1
  99.       IRS=M
  100.       ICS=1
  101.       IMEND=MN-M+1
  102.       MN=M
  103.     6 IOPE=MOD(IOP+3,3)
  104.       IF(IABS(IOP)-3)7,7,2
  105.     7 IF(IOPE-1)8,18,8
  106. C
  107. C        INITIALIZE SOLUTION OF TRANSPOSE(T)*X = A
  108.     8 MEND=1
  109.       LLD=IRS
  110.       MSTA=1
  111.       MDEL=1
  112.       MX=1
  113.       LD=1
  114.       LX=0
  115. C
  116. C        TEST FOR NONZERO DIAGONAL TERM IN T
  117.     9 IF(T(MSTA))10,3,10
  118.    10 DO 11 I=MEND,MN,ICS
  119.    11 A(I)=A(I)/T(MSTA)
  120. C
  121. C        IS M EQUAL 1
  122.       IF(MM1)2,15,12
  123.    12 DO 14 J=1,MM1
  124.       MSTA=MSTA+MDEL
  125.       MDEL=MDEL+MX
  126.       DO 14 I=MEND,MN,ICS
  127.       DSUM=0.D0
  128.       L=MSTA
  129.       LDX=LD
  130.       LL=I
  131.       DO 13 K=1,J
  132.       DSUM=DSUM-T(L)*A(LL)
  133.       LL=LL+LLD
  134.       L=L+LDX
  135.    13 LDX=LDX+LX
  136.       IF(T(L))14,3,14
  137.    14 A(LL)=(DSUM+A(LL))/T(L)
  138. C
  139. C        TEST END OF OPERATION
  140.    15 IF(IER)16,17,16
  141.    16 IER=0
  142.       RETURN
  143.    17 IF(IOPE)18,18,16
  144. C
  145. C        INITIALIZE SOLUTION OF T*X = A
  146.    18 IER=1
  147.       MEND=IMEND
  148.       MN=M*N
  149.       LLD=-IRS
  150.       MSTA=MM
  151.       MDEL=-1
  152.       MX=0
  153.       LD=-MM1
  154.       LX=1
  155.       GOTO 9
  156.       END
  157.